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Abstract 

We present a swift walk-through of our recent work that uses machine learning to fit 
interatomic potentials based on quantum mechanical data. We describe our Gaussian 
Approximation Potentials (GAP) framework, discuss a variety of descriptors, how to 
train the model on total energies and derivatives and the simultaneous use of multiple 
models of different complexity. We also show a small example using QUIP, the software 
sandbox implementation of GAP that is available for non-commercial use. 
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INTRODUCTION 


Molecular scale simulation is a mature field with a wide range of electronic structure methods 
that approximate the solution of the Schrodinger equation in a systematic fashion. For 
larger scale computations empirical interatomic potentials are used, which are nowadays fit 
to data generated by electronic structure models. Together these play a significant role in 
understanding processes on the microscopic level, complementing experiment and theory. 
Computer simulations are regularly used to interpret experimental results and to predict 
properties of materials. 

The power of atomistic simulations would be enormously enhanced if the interatomic 
potentials used to simulate materials were not limited by the simple empirical functional 
forms but accurately approached the Born-Oppenhcimer potential energy surface, similarly 
to the case of small molecules for which quantum chemists have been fitting accurate potential 
energy surfaces for decades. The challenge in the materials field is that rather than fitting 
the total energy of a fixed number of atoms, the task is to find a unique local functional that 
describes the energy of a single atom or bond given its neighbour environment. This local 
energy function must naturally allow for bond forming and bond breaking, i.e. the change 
in the number and identity of the atoms comprising the neighbour environment. 

A number of groups—many of them contributing to the present volume—have started 
research programmes to address this problem using advances in the synthetic understanding 
that recently emerged in statistics and machine learning. These fast-growing fields are con¬ 
cerned with classification, regression and probability density estimation on large and noisy 
data sets, and also with finding suitable variable transformations that allow increased perfor¬ 
mance in these tasks. There are a number of closely related computational frameworks that 
are widely used, including artificial neural networks, stochastic processes (e.g. Gaussian pro¬ 
cesses) and regularised non-parametric optimisation. In this tutorial introduction we focus 
on a particular exposition that allows a succinct presentation of the formalism and how it 
can be brought to bear on the problem of fitting potential energy surfaces for materials based 
on data computed by electronic structure methods. For detailed derivations of the necessary 
fundamental results we refer the reader to the machine learning and statistics literature 1 2 . 
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neighbourhood 

Set of nearby atoms whose positions constitute the input to the local 

energy function evaluated for a given atom. 

descriptors 

Transformation of the positions of atoms in the neighbourhood, obeying 

the desired symmetries of the energy function. Also called features. 

kernel 

Similarity measure between two neighbourhoods, equivalent to the co- 

variance of the corresponding two local energy values. 


Table 1: Definition of central concepts used in fitting accurate potentials for materials. 

METHODOLOGY 

The hallmark of an interatomic potential is that the total energy, E, of a set of atoms is 
written as a sum of range-separated terms, 

E = EE ef + long range contributions (1) 

a iGa 

where ef are local energy functionals with compact support within a radius r cut , and by 
“long range contributions” we mean electrostatics including polarisability, van der Waals 
interactions etc. This is an uncontrolled approximation, since there is nothing about the 
Schrodinger equation that tells us a priori that its solutions can be written in this form: 
the level of accuracy and its applicability in any particular situation has to be tested by 
numerical experiments. The index a denotes the type of contribution: the arguments of a 
local energy term may be any suitable descriptors , e.g. atom-pair distances, bond angles, 
or indeed the complete atomic environment, and the index i counts the instances of these 
terms in a particular configuration, e.g. all bonds for a pair term, all angles for a three-body 
angle-dependent term, or all atoms for an atom-centered term. We can think of descriptors 
as functions that transform the Cartesian coordinates of the atoms in the neighbourhood of 
a given atom. 

In this paper we will only discuss the local energy contribution, although it is clear that for 
many materials in which atoms acquire significant partial charges or have easily polarisablc 
electrons it must be complemented by long range terms. These can either remain completely 
empirical, but may also include parameters that are fitted to data using approaches similar 
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to what are used for the local term. 


Gaussian Process Regression 

We first consider the case of a single type of local energy functional. Using a set of arbitrary 
basis functions {(fh}h=i that take as their arguments any descriptor dj of the neighbour 
environment of atom i, we write the atomic energy e* as 

£i = e(dj, w) = w h cj) h (di), (2) 

h 

where w is a vector of weights Wh corresponding to the basis functions, to be determined 
by the fit. If the prior probability distribution of the weights is chosen to be Gaussian with 
zero mean, i.e. P( w) = Normal(w; 0, <7,„I), the covariance of two atomic energies is 

(£i£j) = ( ^WfcWfc'0h(di)0 h #(dj) j = y^(w /t wy)0/ t (d < )0v(dj) = a 2 w 0 /l (d i )0 /l (d j ) (3) 

\ hh' / hh’ h 

where we exploited that ( WhWh>) = Shh'&w- The inner product of the basis functions in the 
last expression defines the kernel or covariance function 


(4) 

h 

Kernel functions in this application are to be understood as similarity measures between 
two atomic neighbour environments. Every basis set induces a corresponding kernel, and as 
seen below, only the kernel is required for regression, we never need to construct a basis set 
in the space of descriptors explicitly. General requirements on kernel functions are in the 
literature 1 3 . 

Our goal is to predict the energy of an arbitrary atomic configuration, based upon a 
data set of previous calculations. For any set of microscopic observations t —which could 
be the local atomic energies or the total energies of all atoms in a set of configurations— 
the covariance matrix is defined as C = (tt T ), and its elements can be computed using the 
previously defined covariance function. The prior probability of observing t is also Gaussian, 


P( t) = Normal(t; 0, C) oc exp 



(5) 
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The predicted value y of a new test configuration, given previous observations t, has the 
probability distribution 

= (0 

which is also Gaussian. We take the mean of this distribution as the prediction, which can 
be expressed ^ as 

y = k T C“H (7) 

where k is the covariance vector of function values: k = (y t). This shows the real power 
of the Gaussian process approach: the original basis functions we started with and their 
corresponding unknown weights are never required explicitly, the predictions only depend 
on the kernel function C and the previous observations t. 

It can be shown 4 ! that a two-layer neural network with infinite number of hidden nodes 
and hyperbolic tangent switching function is equivalent to a Gaussian process with 

C(&i, cl,-) oaV— |dj — dj| 2 . (8) 

Extra layers in neural networks with more than two layers can be regarded as performing 
a nonlinear transformation on the input coordinates, before the output layers carry out the 
regression task. 

Yet another equivalent approach for fitting functions is kernel ridge regression, where the 
unknown function is expanded as a linear combination of radial basis functions, 

/(d) = £>C(d,d f ), (9) 

i 

and the weights a are optimised by minimising the cost function 

i = E( f ‘-/( d ')) 2 + A ii«n 2 - ( 1Q ) 

i 

If we define the norm as 

| |a| | 2 = «C ck, (11) 

kernel regression is also equivalent to a Gaussian process. The here kernel has the dual role 
of defining both the basis functions and the norm of the weights in the loss function. 
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Total energies 


Atomic energies are unavailable in quantum mechanical calculations, which only provide the 
total energy and its derivatives. From these, we have to predict the local energies. It is 
straightforward to modify equation ([3]) to express the covariance of the total energies of two 
set of atoms, N and M, 

{E n E m ) = ^^Wh'M^Mdj)\ = 

\i£N j£M / \i£N j&M hh' / 

^(whWh'jMdi^h'idj) = <J 2 W M d i)M d j) = (T l^2Yl c '( d *’ d i) 

i&N jeM hh 1 i£N j£M h i£N j&M 

( 12 ) 

Derivatives 

The total quantum mechanical energy of a configuration depends on the relative positions of 
the atoms and, in case of condensed systems, also the lattice parameters. Denoting a general 
coordinate by £, the partial derivative of the total energy is related to the force as 


or to the viral stress as 


, dE OE 

Jka = -w- = -yvr h 4 = r ka 

dr ka dt; 


dE dE , 

Val> = aKr^^ aP 


(13) 


(14) 


where ry Q is the a-th component of the Cartesian coordinates of atom k and h a p is an 


element of the deformation matrix H of the lattice vectors. Differentiating equation (12) 
with respect to an arbitrary coordinate h of configuration N results in 

ddj 

d£h 


dE N \ d(E N E M ) 


' T -EE v * c '( d ” d i 

i£N j£M 


(15) 


S& / 

If is the x, y , or z component of the position of atom k, becomes exactly zero if 
the pair distance |— r*| is beyond the cutoff of the environment, so the hrst sum need 
not be done over all atoms in the configuration. Similarly, the covariance of two derivative 
quantities may be written as 

'dE N dE M \ d 2 (E N E M ) 


dik dxi 


dhdxi 


= <J n 


££^(v*c'(d i ,d J )v;/ d ’ 


i£N j&M 


dh 


dxi ’ 


(16) 
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where the elements of the Jacobian are 


(V d ,C(d„d j )Vl)o (J = 


3 2 C(di,d,) 


(17) 


ddiaddjp 

The local energy e is still predicted by using equation ([7]), but the elements of y are total 
energies or derivative quantities, and the elements of the covariance matrix C are therefore 


computed by equations (12), (15) or (16). The elements of k are the covariance between the 
local energy that we wish to predict and the data that we have available, (e E) or (e dE /d$) 
as appropriate. 


Multiple models 

Interactions in some atomistic systems might be partitioned using a many-body type ex¬ 
pansion - indeed, many traditional interatomic potentials are based on a few low-order 
contributions, such as two- and three-body energies®. We now describe how such models 
can be fitted using Gaussian process regression. For example, truncating the local part of 
equation ([l]) at three-body contributions, the total energy is approximated as 

E = E f + E J 3) < 18 ) 

p G pairs t G triplets 

where and e® are general two- and three-body energy functions, respectively, and pairs 
and triplets in this context may refer to atoms as well as entire molecules. Two independent 
Gaussian processes are used, 


£ (2) (-) = £-i 2 Vi 2) (") 

h 

(19) 

£ (3) (.'.) = (.'.)> 
h 

(20) 


(21) 


where •• anddenotes generic geometric descriptors of pairs and triplets (in case of molecules, 
the descriptors need to describe the whole dimer and trimer configuration). The prior distri¬ 
butions of the two weight vectors are independent Gaussians, so the covariance of the total 
energy of two configurations N and M may be written as 

(E N E M )=a 2 wK Y W C [2) (j>,q) + o 2 wm £ £ C W (M), (22) 

pGpairs jy ^Gpairs^ £Gtriplets^ wGtriplets M 


7 






where we applied the same kernel trick as above and exploited that (w^wffi) = 0 for any 
h and h'. As the two-body terms can, in principle, be included in the three-body terms, 
splitting them appropriately might require setting the the variances of the two terms carefully. 
For example, if 80% of the total interaction energy is expected to be due to pair interactions, 
this information can be built into the prior by using the ratio cr w ( 2 ) : <J W (3> =4:1. 


Compact support 

The local energy terms need to have compact support to be computationally efficient, and 
this is typically achieved by using an explicit spatial cutoff function. In machine learning 
models for materials, the cutoff may be built into descriptors®, so only neighbours within a 
predefined radial distance of the central atom are considered. Alternatively, cutoffs may be 
implemented in the kernels. Consider the pair energy model 

£ (2, (") = /c ».(••) 2>£M 2) (-) (23) 

h 

where / cut is defined such that it goes smoothly to zero as a function of the geometric 
attributes of the pair (e.g. as the distance between them approaches a limit, in case of a 
pair of atoms). The resulting covariance function is 

(e (2) (p)£ (2) (g)) = or 2 w(2) C {2) (p,q)f cut (p)f cut (q). (24) 


In our implementation we use 

' 

1 for r < r cut — d 


/cut(r) = < 


[cos (7T r " r 7+ d ) + 1] /2 


for r cut - d <r < r cut 


0 for r > r cut 


as the cutoff function, where d is a parameter that determines the width of the cutoff region, 
typically 1 A. 


Data noise 

A configuration’s total quantum mechanical energy and its derivatives can be extrapolated to 
exact numerical values, provided the various convergence parameters of the applied quantum 





mechanical method are used appropriately. However, we should still regard these as noisy 
observations when trying to fit a model, for the following reasons: 


(i) the separation into a sum of local contributions is an approximation, 


(ii) our model is additive over various contributions with unknown individual ratios, 


(iii) our model employs a finite cutoff, 


(iv) the quantum mechanical calculations may not be fully converged^ 


Thus we modify our model in equation (JTj) to include a Gaussian noise u E , P(u E ) = 

Norma^z/^; 0, a E ) 

£ = EE £ ? + I * < 25 ) 

Oi i 

and similarly, derivative quantities are modelled as 


dE 

dik 


dik 




(26) 


where P(z^) = Normal(z/g; 0, cr^). As a consequence, equation (12) is modified to give the 
covariance of observed total energies 


(E n E m ) = o 2 w ^ ^ C(di, d j) + ct 2 e Snm, 

i£N j£M 

and the covariance of observed derivatives becomes 


d 2 (E N E M ) 

dikdxi 


<9d7 


= cr„ 


ddj 


E E ^(Vd,C(d i . d j )Vi J )|=l + a(S NM 5 ikX , 

ieN j&M ^ k X - 1 


Sparsification 


(27) 


(28) 


It is easy to see that computing covariance matrices and vectors can become quite expensive, 
especially if derivative quantities are also included. This, combined with the assumption 
that atomic neighbourhood environments are often repetitious, leads to the idea that sparse 
Gaussian processes might be applied. Sparsity is a central concept in machine learning, 
and sparse Gaussian processes are described in detail by, for example, Quinonero-Candela 

*This could be an advantage, as we do not need fully converged quantum mechanical data. 
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and Rasmussen/ or Snelson and Ghahramani 8 . In our adaptation of sparsification, we use 
representative atomic neighbourhood environments, or pairs and triplets etc. The model is 
built using all observations in the dataset and it can be regarded as a projection onto a 
subset of data points, the sparse representation. 

Let us consider a set of configurations, each of which contains an arbitrary number of 
atoms and the corresponding set of total energies, derivatives or both. The observables are 
collected in the vector t. We select a set of environments, the sparse set S, and compute the 
covariance matrices: 

(C 5 s) ss ' = {£ s £ s '), where s, s' G S, (29) 

(Csr)st = {^sEt)-, (30) 

where s G S and t is an index of total energies in t, and 

(C ST ) sr = <^>, (31) 

where s G S and r denotes derivative observables in t. The predicted value at an arbitrary 
atomic neighbourhood environment d* can be calculated from 

£*(d*) = kJ(C55 + Cst^-ttCts) 1 C5 r A T yt, (32) 

where (k*) s = (e*e s ), and A TT is a diagonal matrix, where each diagonal element is o% or 

cr|, depending on the type of observable. As configurations may contain different numbers 
of atoms, we scale a 2 E accordingly. Note that the part multiplying k* from the right is 
precomputed at the training stage, so only k* needs to be computed for each prediction, 
and this scales linearly with the number of sparse points (and not with the total number 
of original data points!). Derivatives of e* are readily available analytically, using the ap¬ 
propriate covariance functions. In practice, we found that the sparse covariance matrix Css 
should be regularised by adding a small positive constant crater to the diagonal values. The 
numerical value of the constant should be as small as possible, without compromising the 
positive definiteness of Css■ Normally crater is 6-9 orders of magnitude less than the diagonal 
elements. 
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Descriptors 

The success of applying machine learning techniques to fit potential energy surfaces depends 
to a large extent on representing the atomic environments appropriately. Transformations 
of atomic positions to which the local energy is invariant, i.e. rotation and inversion of an 
environment about its centre, and permutation of identical atoms should be explicitly built 
in. We presented a detailed study on representing chemical environments elsewhere 6 in which 
focussed on atom-centred neighbourhood environments. Here we describe a few other types 
of descriptors. 

Pairs and triplets 

Pairs of atoms are simply described by the distance between them, but in case of triplets 
the distances need to be symmetrised. If atoms j and k form a triplet with i as the central 
atom, a possible descriptor can be the vector 

Vik + Tij, {r ik - rij) 2 , r jk \ . (33) 

As we mentioned earlier, the covariance function must be augmented by a cutoff function. 
We use /cut ) for the pair terms, and in case of triplets we use / cu t(^ij)/cut(^fc)- 

Water dimers 

It is clear that our approach to symmetrise distances in case of three-body descriptors will be 
overly complicated if we attempt to apply it on more than a couple of atoms. For example, 
the potential energy surface of water can be modelled very accurately using a many-body 
expansion of interactions between water molecules] 9 10 The two-body term in the expansion 
necessitates a descriptor for the water-water dimer, for which we used the pairwise distances 
between the constituent atoms, 15 in total. However, this descriptor in this form is not 
invariant to permuting atoms of the same element. If exchange of hydrogen atoms between 
different molecules is not permitted, the following permutations P that operate on the order 
of atoms must be taken in account: 

(i) swaps of water molecules in the dimer (2) 
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(ii) exchange of hydrogen atoms within each molecule (2x2), 


so 8 in total. Instead of modifying the descriptor, we enforced permutation symmetry at 
the level of the kernel function. If we take an arbitrary kernel function, C(d, d'), that takes 
vector arguments, we can generate a permutational invariant kernel as 


C(d,d') = £c(d,Pd'), 


(34) 


which must be normalised® 


C"( d, d') = 


C"(d,d') 


V , C'(d,d)v , C'(d',d')' 

We used the squared exponential as our starting kernel in the case of water molecules. 


(35) 


SOAP 

We note that our previously introduced® kernel based on “Smooth Overlap of Atomic Posi¬ 
tions” may be interpreted from the function-space view we used throughout this manuscript. 
We represent the atomic neighbourhood of atom i by the neighbourhood density function 
(for illustration, see Figure [I]) 

neigh. 

d*( r ) = exp 

3 

and £j, the atomic energy of atom i can then be regarded as a functional of pi 

£i = e[pi} = I w(r)pi(r)dr (37) 

where the prior distribution of the weights is Gaussian, so 

(w(r)w(r')) = 5(r - r ')<j 2 w , (38) 

resulting in the covariance of two atomic energies 

t i ’(r)pi(r) , u;(r / )pj(r / ) dr dr'\ = cr^ J p^pjir) dr. (39) 

It is useful to note here that if C is a valid kernel, then \C\ P is also valid 3 . This covariance 
function \C\ P is not invariant to rotations, but we may convert it in a similar fashion to what 
we did in the case of the water dimer: 

C\pi, Pj ) = J\C( Pi ,R Pj )\*dR, (40) 


C{pi,pj ) = (£i£j) = / 


i 


2rr 2 

^ u atom 


(36) 
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which must then be normalised, so the final result is 


C"( P i, Pj ) = 




V C '(Pi’Pi)V C '(pj,Pj) 


( 41 ) 


In practice, we evaluate the SOAP kernel numerically by first expanding equation (36) 
in a basis set 

Pi ( r ) = XI c nim9n(r)Y lm ( r), (42) 


nlm 


where Cyj m are the expansion coefficients corresponding to atom i, {g n { r )} is an arbitrary 
set of orthonormal radial basis functions, and Y) m (f) are the spherical harmonics. We form 
descriptors from the coefficients by computing the power spectrum elements 


w = Vc w (c {i) r 
rnn'l / j u 'nlmV-'n'lm) ’ 

m 

and the rotationally invariant covariance of atoms i and j is given by 

C\pi,pj) = Pnn'lPnl'l’ 

n,n' ,1 


(43) 


(44) 


which we normalise according to equation (41). The normalisation step is equivalent to 


normalising the vector elements so C" is, in fact, a dot-product kernel of vectors 

p(0/|pb)j anc [ p (i)/|pO')|. Note that it is often useful to raise C" to a power £ > 1, in order 
to sharpen the difference between atomic environments. To see the details of the above 
results, we refer the reader to our earlier work. 6 n . 


SOFTWARE 

The ideas presented in the Methodology section are implemented in the QUIP package, 
which can be downloaded from the git repository at https://github.com/libAtoms/QUIP. 
Code related to GAP prediction can be obtained under an non-commercial licence from 
http://www.libatoms.org/gap/gap_download.html. Users who wish to use the training 
code should contact the corresponding author. 

QUIP is a molecular simulation sandbox written in object-oriented FORTRAN95/2003, 
with interfaces to python (compatible with ASE), and various other simulation packages, 
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such as LAMMPS, CP2K, CASTEP, and others. QUIP is essentially a collection of objects 
and interfaces that contain and manipulate atomic configurations and interatomic potentials. 

The GAP implementation provides over 20 different descriptors, which can be used with 
two types of covariance functions, the squared exponential 

C(d h dj) = 5 2 exp ( fi?Q g Q djQ ) j ( 45 ) 

and the polynomial kernel 

C'(d i ,d i ) = <5 2 (di-d j + ^) c . (46) 

To demonstrate the training, we provide a simple example, where we train GAP to 
reproduce the well-known Stillinger-Weber (SW) potential 5 for silicon. We used a database 
of Si n silicon clusters, n — 7... 13, sampled from a 2000 K molecular dynamics simulation. 
We used 600 configurations in total, with total energies and forces calculated using the SW 
potential. We used a combination of two- and three-body interactions, with a cutoff of 4.1 A. 
We used the command line 

teach_sparse at_file=data_Si_SW.xyz descriptor_str=-[ \ 

distance_2b cutoff=4.1 n_sparseX=250 covariance_type=ard_se theta_fac=0.5: \ 
angle_3b cutoff=4.1 n_sparseX=500 covariance_type=ard_se theta_fac=0.5} \ 
default_sigma={0.001 0.05 0.01} sparse_jitter=l.Oe-8 e0=0.0 

Table [2] summarises the command line arguments used in this example. The command line 


at_file 

contains the database configurations, in concatenated XYZ hies 

descriptor_str 

parameters of descriptor(s) 

default_sigraa 

the assumed standard deviation of the errors, {cr energy <7f orce cr viral } 

sparse.jitter 

regularisation of the sparse covariance matrix, Opter 

eO 

baseline of atomic energies 


Table 2: Overview of the basic command line arguments of teach_sparse. 


argument descriptor_str contains the parameters of the descriptor, which depend on the 
type. We define the number of sparse points n_sparseX and the type of covariance function. 
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theta_fac is the simplest way to control the 9 a length-scale parameters in the squared 
exponential covariance function: the range of the descriptor values in each dimension a is 
scaled by the constant theta_fac. More descriptors can be concatenated, separated by 
the : symbol, resulting in the fitting of a model that is the sum terms each based on one 
descriptor. 

It is possible to specify an existing QUIP potential as a baseline, so energies/forces/virials 
are subtracted from the target values before fitting. These are added back automatically 
when the potential is called. Naturally, the baseline can be another GAP, resulting in 
hierarchical models with arbitrary level of recursion. 

In the above example, it is possible to check whether the GAP model was able to recover 
both terms of the target model, as they are available analytically. We emphasise that the 
fitting uses total energies and forces only, so the machine learning algorithm has to infer the 
separate two- and three-body terms from this convoluted data. To show the quality of the 
fit, we plot the pair potential in Figure [2] and the angle term in Figure [3] for both the original 
SW and the fitted model. The agreement is rather good in both cases, except at the edges 
of the range where there was no input data. 

The potential hie generated by teach_sparse can be used to compute the total energies 
and similar quantities of arbitrary configurations. For example 

eval at_file=data_Si_SW.xyz param_file=gp.xml init_args={IP GAP} e f 

will compute the total energies (“e”) of the configurations stored in data_Si_SW.xyz as well 
as the atomic forces (“f”) using the fitted GAP model from the hie gp.xml. Another use of 
eval is to compute and print the descriptor vectors for any descriptor type implemented in 
QUIP. For example, 


eval at_file=data_Si_SW.xyz descriptor_str=-[angle_3b cutoff=4.1} 


prints all three-body descriptors, in this case the descriptors defined in equation (33). 
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Figure 1: A two-dimensional illustration of the atomic neighbour density function used in 
SOAP. The white circle represents the radial cutoff distance. 
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bond length (A) 


Figure 2: The two-body term in the Stillinger-Weber potential for silicon. The dashed line 
represents the analytical answer, the solid line shows the GAP fit. The histogram above 
corresponds to the input data to the fit. 
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Figure 3: The three-body term in the Stillinger-Weber potential for silicon. We fixed the 
two neighbours at 2.5 A and 2.8 A and varied the bond angle and plotted the sum of all three 
interactions between the three atoms. The dashed line represents the analytical answer, the 
solid line shows the GAP fit. The histogram above represents the input data to the fit. 
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